MULTIGRID PRECONDITIONING OF LINEAR SYSTEMS 
FOR INTERIOR POINT METHODS APPLIED TO A CLASS OF 
BOX-CONSTRAINED OPTIMAL CONTROL PROBLEMS 
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Abstract. In this article we construct and analyze multigrid preconditioners for discretizations 
of operators of the form "Dx + K*JC, where Dx is the multiplication with a relatively smooth func- 
tion A > and /C is a compact linear operator. These systems arise when applying interior point 
methods to the minimization problem miun i(||ASji — /p + /3||up) with box-constraints u ^ u ^ u 
on the controls. The presented preconditioning technique is closely related to the one developed 
by Draganescu and Dupont in 1131 for the associated unconstrained problem, and is intended for 
large-scale problems. As in |13l . the quality of the resulting preconditioners is shown to increase as 
h X 0, but decreases as the smoothness of A declines. We test this algorithm first on a Tikhonov- 
regularized backward parabolic equation with box-constraints on the control, and then on a standard 
elliptic-constrained optimization problem. In both cases it is shown that the number of linear iter- 
ations per optimization step, as well as the total number of fine-scale matrix-vector multiplications 
is decreasing with increasing resolution, thus showing the method to be potentially very efficient for 
truly large-scale problems. 
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1. Introduction. In this work we present a multigrid preconditioning technique 
for solving linear systems arising when applying interior point methods to the control- 
constrained optimal control problem 

minimize Jp{u) ^\\ICu - ff + ^\\uf, u e U^d , (1.1) 

where 11 C M'^ is a bounded domain, /3 > 0, /C : L'^{il) — > L'^{n) is a linear, compact 
operator, and the set of admissible solutions is given by 

^ad = {u G L^{i}) : u^u^u a.e.}, 

with u, u e i^(il), and u^u. The following examples form the main motivation for 
the present work: 

Example A: Box-constrained time-reversal for a parabolic equation. Con- 
sider the linear parabolic initial value problem 

dfU + Ay — , onl7x(0,oo), 

y{x,t)^0 , on ai7 X (0, cx)) , (1.2) 
y{x, 0) = u{x) , for a; € 17 , 

where ^ is a linear elliptic operator, and denote the solution map by S{t)u y(-, t). 
To formulate a control problem we define /C = S{T), where T > is a fixed "end- 
time" T > 0; the resulting optimization problem is controlled by the initial value u. 
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Note that if the box constraints in are left out, i.e., Z-/ad — L'^i^), then (jl.ip is 
the Tikhonov-regularized formulation of the ill-posed problem 

ICu = f. (1.3) 

Example B: Elliptic-constrained distributed optimal control problem. 

In this example we let 

/C = A^i , 

where A is the Laplace operator acting on HQ{il). In this case the problem (II. ip is 
usually formulated as the PDE-constrained optimization problem (e.g., see Borzi and 
Schulz [1]) 

minimize i \\y - ff + | \\uf 

subj. to Ay — u , y £ H}^{Vl), u^u^u a.e. in il. 

Our primary motivation is rooted in solving large-scale inverse problems like the 
one in Example A, which is a simplified version of the problem considered by Akgelik 
et al. in [2]. There the question was to identify the initial concentration u = y(-, 0) of 
a contaminant released in the atmosphere in a given geographic area (the Los Angeles 
Basin) given later-time measurements at various fixed locations. The spatio-temporal 
evolution of the concentration y of the contaminant is assumed, like in Example A, 
to be governed by an advection-diffusion equation with known wind-velocities and 
contaminant diffusivity that can be formulated as (|1.2p . In Example A we consider 
the case where the measurements are taken at all points in space but at a single mo- 
ment in time T > 0, therefore the data / is the entire state at time T. The problem 
thus becomes to invert a compact operator which, as is well known, is not contin- 
uously invertible. As a result, a naive approach to inverting /C is unstable, in that 
small perturbations in the "measurements" / result in exponentially large errors in 
the solution. If the "exact" measurements / are resulted from applying /C to a "true" 
initial value u, that is, / = /Cm, then various regularization techniques [15) are em- 
ployed so that the computed solution us of the (5-perturbed problem ICug = fs (where 
1/ ^ /(sll ^ <5) converges to the "true" solution u as 6 I Q, the most commonly used 
being the Tikhonov-regularization. One issue not resolved by the Tikhonov regular- 
ization is that the solution us of the regularized problem may exhibit non-physical, or 
otherwise qualitatively incorrect behavior: for example, if the concentration needs to 
have values in [0, 1], it is well known that ug may exceed these limits. In addition, if 
the true initial contamination event is localized, then ug oscillates around zero, and is 
not localized. However, if explicit constraints u = Q,u=l are set, then the solution 
is physically relevant, and is also localized. In Figure [TTT] we show the solution of the 
inverse problem where the target solution (true initial concentration) is represented 
by two localized sources of contamination, with one dominating the other in size. 
The solution of the unconstrained problem (with 5 down to roundoff error and /? opti- 
mized) does not show two sources of pollution, unlike the solution of the constrained 
problem, from which one can clearly identify two separate components in the sup- 
port of the initial value. In addition, in the presence of nonlinearities, with reaction 
terms that are sensitive to the sign of their arguments, the importance of physically 
meaningful box constraints cannot be overstated. A similar situation is encountered 
in image-deblurring, where the target (gray-scale) solution takes values in [0,1]; if 
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Fig. 1.1. Solution of the constrained vs. unconstrained inverse problem, with the target solution 
{initial concentration) being a smooth function supported on two disjoint intervals. 



these boundaries are strictly enforced, then the quahty of inversion (deblurring) is 
significantly improved |34j . 

The unconstrained optimization problem, where Ua,^ = L'^{fl) in p.ip . ultimately 
reduces to solving the normal equations 

{PI + JC*K.)u = K.*f . (1.5) 

For compact operators these systems either represent or resemble very well inte- 
gral equations of the second kind. Starting with the works of Hackbusch [5D] (see 
also [12]) niuch effort has been devoted to developing efficient multigrid methods for 
solving (dH]), e.g., see [HI [33 HI EZl O [13] and the references therein. We should 
point out that multigrid methods were originally developed for solving elliptic equa- 
tions |2TJ [To] , and later significant efforts were devoted to extending these methods 
to other important differential equations such as advection-diffusion and the Navier- 
Stokes equations (e.g., see pjj and the references therein). For elliptic equations, 
typically the goal is to reduce the condition number of the discrete system from 
0(/i"2) to 0(1), which results in a solution process that solves the equation in a num- 
ber of iterations that is mesh-independent. However, the discrete version of (jl.Sp has 
a condition number which is 0(/3^^), with the bound being independent of h. More- 
over, even for /3 = 0, conjugate gradient (used as a regularizer) already solves (|1.5p 
in a mesh- independent number of iterations [f^; in other words, mesh-independence 
is nearly effortless for integral equations. Instead, for systems like (|1.5p multigrid 
is used to further reduce the condition number of the preconditioned system. For 
example, in |13| it is shown that by using specially designed multigrid precondition- 
ers, the preconditioned version of (|1.5p has a condition number bounded by 0(/i^//?), 
the consequence of which is interesting at least from a theoretical point of view: if 
/3 > is kept fixed (which is normally not the case in practical applications) the 
number of iterations required to solve the problem decreases with /i J, to the 
point where only one iteration would be enough to solve the problem with sufficient 
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accuracy. This fact, already known to Hackbusch i20. , constitutes a departure from 
the usual multigrid framework, where, as stated before, the goal is to achieve mesh- 
independence (bounded number of iterations as h ^0). Hence, as a result of specific 
multigrid preconditioning, the solution process for the unconstrained problem (jl.Sp 
requires fewer and fewer fine-scale matrix-vector multiplications (mat-vecs) as /i J, 0. 
The main contribution of the present work is to show that such performance can also 
be achieved in the presence of explicit box constraints on the control, as formulated 
in (fTTl) . 

Multigrid methods have long been associated with solving large-scale problems, 
and beginning with the work of Hackbusch [181 [19], and especially over the last decade, 
significant efforts were concentrated on devising efficient multilevel methods for op- 
timal control problems [31 [3 [HJ |40l |29l . A more detailed discussion of the subject 
and many references can be found in the recent review article by Borzi and Schulz [8] . 
While most - not all - of the aforementioned articles discuss unconstrained problems, 
the addition of box constraints pose additional challenges associated with the presence 
of the Lagrange multipliers related to the inequality-constraints, which in general are 
less regular than the solutions. Optimization methods for such problems with bound 
constraints typically fall in one of two categories: active-set type strategies, especially 
semismooth Newton methods (SSNMs) and interior point methods (IPMs). Over the 
last decade both IPMs and SSNMs have consistently attracted the attention of the sci- 
entific community due to their proven efficiency in solving distributed optimal control 
problems with PDE constraints. Both strategies lead to superlinear local conver- 
gence [38l |39l [21] and lend themselves to analysis both in a finite dimensional setting 
and in function space [ST] [41] (see also [23), which is a critical stepping stone to 
proving mesh- independence [5S]. Each of IPMs and SSNMs consists of an outer it- 
erative process that further requires solving one or two linear systems at each outer 
iteration. For large-scale problems the solution of these linear systems often becomes 
the bottleneck of the computation. In terms of their linear algebra needs, IPMs and 
SSNMs exhibit significant differences and require separate treatment. For SSNMs, 
the linear systems involve a subset of unknowns and equations of the Hessian of the 
cost functional, while for IPMs the systems have the same structure as the systems 
arising in the unconstrained problem, but contain additional terms on the diagonal 
which usually are a source of extreme ill-conditioning. The question of efficient multi- 
grid preconditioning of the linear systems arising in the semismooth Newton solution 
process is the subject of current research [11] . 

The focus of the current work addresses the question of efficient multigrid precon- 
ditioning of linear systems arising in the IPM solution process. Multigrid and IPMs 
have been shown to work well together for elliptic variational inequalities (obstacle 
problems) [3] , and for some classes of problems where the Hessian of the cost func- 
tional is elliptic f4^. Instead, the Hessian of the cost functional in (jl.l|) is a compact 
operator and requires a significantly different approach. 

In this article we treat /C and its discretizations as black-box operators, and we 
employ the first-discretize-then-optimize strategy. We apply specific primal-dual 
IPMs to the discrete version of (jl.ip . and we develop multigrid preconditioners for 
the linear systems arising at each IPM iterate. As shown in Section [S] if the discrete 
optimal control problem is formulated appropriately, then the linear systems to be 
solved involve matrices of the form D + K*K, where D is a diagonal matrix, K is 
the discrete representation of /C, and K* is the adjoint of K with respect to a certain 
discrete inner product. When using standard finite elements, the matrices K, K* 
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are dense; consequently, for large-scale problems, they cannot be formed and the 
systems are solved using iterative methods. Since residuals can be computed at the 
equivalent cost of two applications of K, residual computation is expected to be very 
expensive, therefore efficient preconditioners are critical for minimizing the number 
of necessary mat-vecs. Our strategy in this work is to adapt the multigrid techniques 
developed in [13j . where the matrix D had the form /3I. To analyse the resulting 
multigrid preconditioner we interpret D as being the discretization of an operator 
Vx, where (T)xu){x) — X{x)u{x) is the operator representing pointwise multiplication 
of a function u with a smooth functioro A. From a technical perspective, our main 
accomplishment consists of showing that the operators of the type ICDx together with 
their discretizations satisfy a set of smoothing conditions shown in [T3] to be sufficient 
for the multigrid preconditioner to have the desired qualities. 

For simplicity and concreteness we restrict most of our study to the two-dimensional 
case, and we consider a standard finite element discretization for /C using triangular 
elements and continuous piecewise linear functions. As will result from the analy- 
sis, these techniques can be easily generalized to three dimensions and rectangular 
elements, however, the extension to higher degree finite elements is not obvious and 
forms the subject of current research. 

This article is organized as follows: After formally introducing the discrete opti- 
mization problem in Section [2j we discuss the specific linear algebra requirements of 
the interior point methods in Section [31 Section S] is central to this work as it presents 
the analysis of the two-grid preconditioner, the main result being Theorem 14.91 In 
Section [5] we develop a multigrid preconditioner that preserves the qualities of the 
two-grid preconditioner. Further, we apply the methods to Examples A and B in 
Section [5] and show some numerical results. 

2. Notation and discrete problem formulation. Throughout this paper 
we shah denote by (n) , H"" (n) , H^" (n) (with p e [1, oo], m € N) the standard 
Sobolev spaces, while || • || and (•, •) are the L^-norm and inner product, respectively. 
Let H-"^{n) be the dual (with respect to the L^-inner product) of iJ™(0) n H^{fl) 
for TO > 0. If AT is a Banach space then £(A) denotes the space of bounded linear 
operators on X. We regard square n x n matrices as operators in £(R") and we 
write matrices and vectors using bold font. If A is a symmetric positive definite ma- 
trix, we denote by (u, v)^ ~ v-^Au the A-dot product of two vectors u, v, and by 
|u|a = (u, u)^ the A-norm; if A = I we drop the subscript from the inner product 
and norm. The space of to x n matrices is denoted by Mmxn] if to = n we write M„ 
instead of Mmxn- Given some norm || • \\s on a vector space X, and T G £(<-t'), we 
denote by ||r||s the induced operator- norm 

||T||s = sup ||Tm||s . 
uex, ||tt|U=i 

Consequently, if T e £(L^(ri)) then ||T|| (no subscripts) is the operator-norm of T. 
If A" is a Hilbert space and T e £{X) then T* € £{X) denotes the adjoint of T. 

We assume that C is a bounded, polygonal domain and that /C is discretized 
using continuous piecewise linear functions on triangular elements. We consider the 
usual multigrid framework where the operator is discretized at several resolutions. 
Let 7?io be a triangulation of the domain fl, and define Th/2 inductively to be the 



^This can always be accomplished, A can be any smooth (here is sufficient) interpolant of the 
discrete function representing the diagonal of D. 
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Goursat refinement of Th for a\\ h E I with 

I = {ho/r -.1 = 0,1,2,...} , 

where each triangle in T e 7?i is cut along the three lines obtained by joining the 
midpoints of its edges. Note that {Th)h^i is a quasi-uniform triangulation and the 
usual approximations hold [9]. We define 

V/i = {u £ C{Vl) : VT e Th u\t is linear, and wjoji = 0} , 

so that V/i/2 C V/i C Hq{VI). We should note that zero-boundary conditions for 
the controls are consistent with the examples considered, and present a convenient 
framework for the analysis, but alternate boundary conditions can be considered. For 
Example A, setting u to be zero on dU, is natural, while in Example B the boundary 
values of the control do not enter the discrete problem unless higher order cubatures 
are used in the discretization. If Nh = dim(V/i) and P^, . . . , P^^ are the nodes of Th 
that lie in the interior of n let It ■ C{n) — >■ Vh be the standard interpolation operator 

Nh 

Xhiu) = J2^iPt)^1 , 
1=1 

where (p'^,i = l,...,Nfi are the standard nodal basis functions. Given a family of 
positive weight-functions {w^jh^i C Vh we define the mesh-dependent inner products 

{u,v)^ = J2whiPhu{Pt)v{Pt), for u,veVh, 
1=1 

and let \u\h = \J {u, u)f^. In order to satisfy {■,-)h ~ (''') ^ close as possible we 
replace exact integration on each triangle APiP2^'3 with the cubature 

/ /(.)dx«Q(/) = ^^^/(PO. 
This defines the weight functions Wh 

MPh^l E ^'^^^(^)- (2.1) 

Since the cubature Q is exact for linear functions 1351 we have 

{u,v)^ = / Ih{uv), for aU u,v eVh ■ 
Jn 

Moreover, if the grids are quasi-uniform, then h~^Wh are uniformly bounded and 
bounded away from with respect to h d I, therefore by Lemma 6.2.7 in [9] there 
exist positive constants Ci , C2 independent of h such that 

C'lM ^ Mh ^ C2\\u\\, yueVh. (2.2) 

We should point out that the norm-equivalence (|2.2p extends to operator norms for 
operators in £(L^(ri)), which allows us to interchange ||r|| with |||r|||h when needed as 
long as we factor in a mesh- independent constant. 
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We assume that for each /i e / is given a natural discretization JCh G S,{Vh) of /C, 
so that JC,ICh satisfy the Smoothed Approximation Condition (SAC): 

Condition 2.1 (SAC). An operator M together with its discretization Aih is 
said to satisfy the Smoothed Approximation Condition if there exists a constant C{M.) 
depending on Ai.Q^Tho and independent of h so that 

[a] smoothing: 

^ C{M) \\u\\ , Vu e L^in), m = 0, 1, 2 ; (2.3) 

[b] smoothed approximation: 

\\Mu - Mhu\\H^.^n) C{M)h^~"' \\u\\ Vu e Vh, m 0, 1, /i e / . (2.4) 

Given two discrete functions Uii,Uh G V/i representing u,u, respectively, we now 
define the discrete optimization problem 

minimize j/^^u) ""^^ i|||/C.^ - f.fl + ^\Ml u e U^^ , (2.5) 
where fh G Vh represents / and the set of discrete admissible solutions is given by 
^Id = e V, : Uh{P^^) ^ < ^P^^) Vz - 1, . . . , Nh} . 

The formulation (|2.5p needs a few comments. First we remark that the use of the 
discrete norm ||| • ||| h instead of || • || is essential in order for the linear systems to be solved 
at each outer iteration to have a form amenable to efficient multigrid preconditioning. 
In other words, we have chosen a discretization that allows for an efficient solution 
process. Second, if u and —Ti are convex and continuous (e.g., when they are constant), 
then the choice 

Uf^^Ih{u), Uh=Ih{u) (2.6) 

implies U^^ C Z^ad- This construction can be easily generalized to three dimensions 
and/or tensor-product finite elements. 

We obtain a matrix formulation of (j2.5p by representing all vectors and operators 
using the standard nodal basis functions (/s^, i = 1, • • ■ , Nh. More precisely, if we define 
T : M^" ^ Vh by 

T(u) = ^ Uiip^, where u = [wi, . . . , un^V ; 

then the matrix K/j = T^^IChT, regarded as an operator in £(R^''), represents ICh 
with respect to the nodal basis. If W/i is the diagonal matrix with diagonal en- 
tries (w?i(Pj''))i^iscjVfe, and th,Uf^,Uh represent fh,Uh,Uh respectively, then is 
equivalent to 

minimize J^(u) =^ i|K;jU - f,,|^^ -f ^|u|^^ , u u,, , (2.7) 

where the inequality u ^ v between vectors is meant coordinate- wise. When operating 
on a single grid we will omit the subscript h for matrices and vectors. 

Existence and uniqueness of solutions for both and (12. 5p follows from the fact 
that J7^, Jn are uniformly convex andZ^/ad, ^ad ^'"'^ convex sets (e.g., see Theorem 1.43 



8 



A. DRAGANESCU AND C. PETRA 



in Furthermore, of. Lemma 1.12 in the solution u of (|l.ip is characterized 

by the foUowing condition: there exist A, A e L^(f2) so that 

(/3/ + /C*/C)m + A - A = /C*/ , 

u^u, A^O, A(u~u)=0 a.e., (2.8) 
u ^ 11, A ^ 0, A(u — u) = a.e. 

3. Interior point methods and linear systems. In this section we briefly 
discuss the specifics of the interior point method we use for solving the discrete opti- 
mization problem (12.71) . and we describe in detail the linear systems that need to be 
solved at each outer iteration. If we denote 

A /?W + K^WK , (3.1) 

then after a rearrangement of the terms in the objective function Jp and dropping 
constant terms we write (|2.7|) as a regular convex quadratic problem with affine con- 
straints in standard form: 

minimize ^u-^Au — (K-^Wf)-'"u 

' _ (3.2) 

subj to: u ^ u ^ u . 

Since W is a diagonal matrix with positive entries, the matrix A is positive definite, 
therefore the above problem has a unique solution (see [32], p. 320). Let us denote 
the Lagrangian corresponding to the QP (|3.2p by 

L(u,vi,V2) = iu^Au- (K^Wf)^u- v^(u- u) - vf (u- u), 

where Vi , V2 are vectors of non- negative multipliers corresponding to the inequality 
constraints. Then the gradient and Hessian of the Lagrangian are given by 

Viu(u,vi,V2) = Au-K^Wf + V2-V1, V2Luu(u,v) = A . 

Since the the linear independence constraint qualification holds, the unique solution u 
of (|3.2p satisfies the Karush-Kuhn- Tucker (KKT) conditions 

{Au + ^2 - vi = K^Wf 
u^u, vi ^ 0, vi-(u-u) = 0, (3.3) 
U ^ U, V2 ^ 0, V2 • (u - u) = , 

where Vi, V2 are the multipliers, and u-s denotes the component-wise product. More- 
over, since V^Luu is positive definite, the above KKT conditions are also sufficient. 
The primal-dual IPM consists of solving the perturbed KKT system 

Au + V2 - vi = K^Wf , 

u > u, vi > 0, vi • (u - u) = ^e , (3.4) 
U < U, V2 > 0, V2 • (u — u) = ^e , 

whose one-parameter family of solutions (u(/i), vi (/i), V2(/i)) defines the central path. 
As usual, e = [1, 1, . . . , 1]-^ e R^*". Practical IPM algorithms produce solutions that 
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lie sufficiently close to the central path and converge rapidly to (u, vi, V2). An exam- 
ple of such method is Mehrotra's predictor-corrector algorithm. Initially introduced 
for linear programming (31) . the method was successfully adapted to convex QPs 
and emerged in the last fifteen years as the (arguably) most practical and efficient 
algorithm for this class of problems. For this project we used Matlab to implement 
Mehrotra's method for convex QPs from OOQP (see [16] for details). 

To describe the method we first consider the the linear system defining the Newton 
direction {Su, S-vi, 6^2) for 



A^u + (5v2 — ^vi = K'^Wf - Au — V2 
Vi(5u + (U - U)5vi = ^e-vi-(u-u), 
— V2(5u + (U — U)(5v2 = /ie — V2 • (u — u) , 



Vl 



(3.5) 



where U, U, U, Vi, and V2 are diagonal matrices with the diagonal given by the 
vectors u,u, u, vi, and V2, respectively. In Mehrotra's algorithm, given the current 
iterate (u,vi,V2), one first computes the predictor direction (Ju", 5v°, Jv^) as the 
solution of p.Sp with fi — 0. Secondly, the corrector direction {du,Svi,5v2) is the 
solution of a linear system that differs from (13.51) only in the right-hand side, namely: 



A5u + (5v2 — Svi = K^Wf — Au — V2 + vi , 
Vi5u+ (U-U)5vi = o-^e- Vl • (u-u) + (5u" -Jv", 
-V2(5u+ (U- U)5V2 = cr/ie - V2 • (u - u) + (5u° • (5vf , 



(3.6) 



where /i = ((u — u)-'"vi + (u — u)-^V2) /{2Nh), and cr > is a centering parameter 
that is computed accordingly to Mehrotra's heuristic. Therefore, both the predictor 
and the corrector step involve a system - the augmented system - of the form 



(3.7) 



Since the matrices on the second and third block-rows of (13.71) are diagonal, we can 
substitute 6vi = (U - U)"i(r^,i - Vi(5u) and 6^2 = (U- U)~i(r^2 + V2(5u) into the 
first block-row to obtain the reduced system 



A I I 




" Su - 






Vl (U - U) 




(5vi 






- V2 (u - U) . 




. <^V2 . 




- ^V2 - 



[A + (U - U)-iVi + (U - U)-iV2] 5u^r 



(3.8) 



with 



r = r„ 



(U-U)-ir,, -(U-U)-ir2 



We note that the matrix of the reduced system (j3.8p is symmetric positive definite, 
while the matrix of (13. 7p is similar to the symmetric indefinite matrix 



A I 

I -Vj:1(U-U) 
I 



I 



-V2-i(U- 



U) 



If strict complementarity holds for at least one coordinate in the pair (u,vi) (resp. 
(u, V2)), then the diagonal matrix Vj^^(U — U) (resp. V2"^(U — U)) has increasingly 
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small and/or large entries as the interior-point algorithm approaches the solution. It is 
easy to see that the largest eigenvalue of C is larger than any of the diagonal entries of 
VJ~^(U — U) or V^"'^(U — U), while the smallest (in absolute value) of its eigenvalues 
are 0(1). Therefore, the matrix C, and hence the system (|3.7p is ill-conditioned. 
Several equivalent reformulations of the system (13.71) are proposed in [4] (also see [T7]V 
As in [21 [13] , in this article we use the reduced form ()3.8p which, for the problem under 
study, is symmetric positive definite. While p.8|) also suffers from the well known 
ill-conditioning of interior point methods, and in fact is even more ill-conditioned 
than p. 71) if /? ^ 1, the condition number of p.8|) can be significantly reduced first 
by rescaling both the unknowns and the equations (left- and right- preconditioning), 
as shown below, and then by two-grid and multigrid preconditioning, as discussed in 
Sections [H and [5] . 

Given a vector m — [mi,m2, ■ ■ ■]'^ we denote by Dj^i the diagonal matrix with 
entries Da = nii, and by p • /m the vector [pi/'mi,p2/m2, . . .]'^ . With A as in (|3.ip . 
w = w/j = [wi, W2, • . . ,WArJ^, and 

m = vi • /(u — u) -I- V2 • /(u — u) , 

the reduced system p.8p can be written as 

(Dm+p^ + K^WK) du^r. (3.9) 

Left- multiplication with W^^ further yields 

(D(^/w)+^e + W-iR^WK) Su = W-ir . (3.10) 

Let p = y/ (m/w) + /3e (component-wise). By rescaling Su' = DpJu , and factoring 
out Dp in (|3.10|) . the system becomes 

(1 + W-^L^WL) (5u' = Di./pW-^r , (3.11) 

where 

L = L/i = K/i Di./p , 

and we used the commutation of the diagonal matrices W^^ and Di./p. We prefer 
to write p. lip in compact form as 

(I H) 5u' = r' , (3.12) 

with 

H = H„ =■ (W,0-1l5^W„L„ , 

and r' = Di./pW~^r. Note that the matrix in (13. 9p is symmetric positive definite, 
while the matrix H is symmetric (and positive definite) with respect to the W-dot 
product. Furthermore, it is easy to see that ||H|| = 0(/3^^||Kp) which implies that 

cond(I + H) = 0(/3-i||Kp) , 

independently of the mesh parameter h. However, for the model problems considered, 
the matrix (I + H) is dense (in standard representation), and therefore for large-scale 
problems it cannot be formed and/or stored. Matrix-vector multiplication can be 
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performed at a cost equivalent to two applications of the matrix K, hence residual 
computations are expensive. So (|3.12p has to be solved using iterative methods, 
and for increased efficiency we need high-quality, matrix-free preconditioners. As 
it turns out, it is the system p.l2|) rather than (|3.9p that renders itself to good 
multigrid preconditioning. In the next sections we develop a multigrid preconditioner 
for (|3.12p under the assumption that m represents a positive and relatively "smooth" 
function /i/j. 

4. The two-grid preconditioner. In this section we develop and analyze a 
two-grid preconditioner for the linear system p.l2p . The work relies on the multigrid 
techniques developed by Draganescu and Dupont in for (|1.5p . As will be shown, 
the constructed preconditioner has, under certain hypotheses, optimal order quality 
with respect to the discretization parameter h. 

4.1. Algorithm design. For the purpose of algorithm design and analysis it 
is advantageous to regard (|3.12|) as an equation in Vh rather than E^'', so we have 
to identify the operator in £,{Vh) that is represented by H^. First we define for 
A G -L°°(ri) the multiplication- by- A operator 2?a by 

Vxu — Xu , 

and its discrete version V'^ E £(V/i) by 

V^^u = IhVxu . 

Given a vector m G M^'' we define a function /i/i e Vh by setting /i/i(P/') — nii, 
i = 1, . . . , A^h; it follows that the diagonal matrix D(„i/w)+;3 represents the operator 

To simplify notation let 

\h^{^ih/wh) + p . (4.1) 

Then L/j = K/jDi./p represents the operator 

Ch - . (4.2) 

If we denote by C*j^ the dual of Cu with respect to the (•, •)^-inner product, that is, 

then C*f^ is represented by W^"^L^W/j, so H/j = (W^)^^L^W/jL;j represents the 
operator 

Hh ■ (4.3) 

Hence, the operator we need to invert for solving (|3.12p is 

gh = i + nh. (4.4) 

Note that the operator Gh is symmetric with respect to (•, •)^, i.e., Qh — G^- 
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The idea behind the proposed two-grid preconditioner for Qf^ lies in the "smooth- 
ing" properties of More precisely, we regard Ch as a discretization of 

for some function A for which 

If A (assumed to be > /3 > 0) is relatively smooth (e.g., it can always be chosen to 
be C^), then the multiphcation operator T^Yjy/x neither smoothing nor roughening, 
so the fohow-on application of /C results in smoothing. An alternative point of view is 
that T)^!^ is bounded in £(L^), and if /C is compact, then ICD^^^ is also compact. 
Hence, it is natural to assume that "Hh = C-h^h is "smoothing" , even though £^ has 
no direct connection with the dual of JCD^^^ in -C(i^). 

We consider the L^-orthogonal splitting of the discrete space 

Vh = V2h ffi mh , (4.5) 

and let tt = TT2h be the L^-projector onto V2h- Following [33J [231 [121 US] , '^e propose 

Mh = p + G2h-K (4.6) 

as a two-grid preconditioner, where p = p2h ~ I — T^2h is the projector on yV2/t, and 
the coarse function X2h entering the definition of G2h is given by 

X2h = IhXh ■ (4.7) 

The operator Mt can also be regarded as an additive Schwartz preconditioner with 
respect to the splitting (|4.5p (see [23]). Moreover, the inverse of Nh is given by 

5,/= AA-^=p + g->- (4.8) 

For developing a multigrid algorithm of comparable quality with the two-grid pre- 
conditioner we follow the same strategy as in [T^l [13], which we briefly outline in 
Section [5] 

As shown in the analysis, the use of the L^-projector 'K2h in the definition of 
as opposed to other projectors or restriction operators turns out to be critical for the 
quality of the preconditioner. Unfortunately 'K2h is not symmetric with respect to 
(•, •)^, therefore Mh and Sh are symmetric neither with respect to (•, •)^ nor to (•, ■)^2, 
but they are almost symmetric. At the root of this problem lies the use of mesh- 
dependent norms in the formulation of the discrete optimization problem (j2.5l) , which 
in turn was necessary for the linear systems inner to the interior point algorithms to 
be of the form p.9p . that is, to have a diagonal matrix Dm+,3w added to K^WK. 
Had we used the exact L^-norm in the discrete formulation (|2.5p . then the matrix in 
the system p.9p would have had the form Dm -I- /3M + K-^MK, where M is the mass 
matrix, and this form is less convenient for preconditioning. 

In order to describe the matrices representing and Sh we consider the prolon- 
gation operator 3^ € ^Nhxn^h representing the natural embedding of V2/1 into Vh, 
and we define the restriction I{.2h G ^'^N2h-xNh by R.2h = ^^''J^, d being the dimension 
of the ambient space (here d = 2). Then ■K2h is represented by the matrix 

U2h = M.-^, ■ R2h • Mft , 
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where (resp. M2/1) is the rescaled mass matrix on the fine (resp. coarse) mesh, 
defined by (M^)^^ = h^"^ {f^: v'j)- Note that 112^ is a matrix of type x ^/i, and 
that M.2h — R2/1 ■ • J/i. Furthermore, the square Nh x A^^ projection matrix is 
given by Ph = J/i • 'n.2h, so that Pf^ = P^- The projector p (/ — tt) is represented 
by the matrix Q/j = (I — P/J e Mn^xNh- Finally, Sh is represented by 

Sh = Qh + JhG-^U2h ■ (4.9) 

Oftentimes in practice the exact projection Il2h in (|4.9p is replaced by the restric- 
tion R2/i, and Qh is taken to be (I — R2/i). While for certain problems this is a viable 
option , for the applications considered in this work the quality of the preconditioner 
is significanly diminished by this change. 

4.2. Algorithm analysis. Our analysis consists of evaluating the quality of the 
two-grid preconditioner by eventually estimating the spectral distance d„{Sh,G^^) 
(Theorem l4.9p . and is performed over three steps. First we evaluate the norm-distance 
\\Qh — A/ft II under the assumption that Hh satisfies Condition 14.11 below. Second, we 
show that if C and Ch verify Condition 12.11 then Hh satisfies Condition 14.11 Third, 
we show that C,Cii satisfy Condition 12. II fSAC) with a constant C{£) depending on 
A and on C(JC), that is, the constant associated to /C,/Ch satisfying SAC. For specific 
applications, the fact that /C,/Cft satifsy SAC is normally verifiable, as is shown in 
Section [51 

Condition 4.1. The operators Hh are symmetric with respect to (•, •)^, positive 
semidefinite, and uniformly bounded with respect to h Cz I , that is, there exists a 
constant C{'H) > independent on h such that 

tHhWh^CiH) , Vhel . (4.10) 

Moreover, there exists p > so that 

WCHh - n2hTT2hH < C{n)hP\\u\\, for all u e V,,, hel . (4.11) 

For linear splines the optimal approximation order is p = 2, but for certain problems 
and discretizations the actual rate may be suboptimal. 
Lemma 4.2. IfHh satisfies Condition 14.11 then 

\\gh-^fh\\^c{H)hp . (4.12) 

Proof. This result is an immediate consequence of (j4.11l) . since 

Gh - Mh = {I + Hh) -ip+{I + n2h)7T) = Hh - H2h7T . □ 

Verifying that (|4.1ip of Condition 14.11 holds under some general hypotheses is non- 
trivial for this problem due to the presence of multiple inner products that have to 
be taken into consideration. More precisely, if C*^ were the dual of £h with respect 
to (•, •) instead of (•, then Condition 14.11 would follow from the approximability of 
C by Ch together with the smoothing properties of £, as is shown in |13j (proof of 
Theorem 4.1). Hence a natural requirement is that {■,-)f^ approximates (•, •) well. 

Lemma 4.3. With Wh chosen as in (j2.ip there exists a constant C = C{Tq) > 
independent of h such that 

\{u,v)h~ ^ Ch^\\u\\Hi(n) ■ \\v\\m{n) , "iu.v £Vh- (4.13) 
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Proof. By Theorem 4.4.4 in P], given d — 2, there exists a constant C > that 
depends on the "chunkiness" of the initial triangulation To but is independent of h 
and oi T G Th, so that for aU u,v E Vh 

II '-7- / \ II ^ ^7 2i I linear on T o I t—j i 

After summing over ah triangles T G T/i we obtain 

C/i^||Vm • Vu||li(o) ^ C/i^|'u|^i(f^) • |w|Hi(n) ■ 

The conclusion follows from 



\{u,v) - (m,w),J = 



{uv - Ih{uv)) 



^ \uv -Ih{uv)\LH^n) . □ 



Throughout this section we denote by C{C) a generic constant that is propor- 
tional to the constant C{C) of Condition 12.11 if the proportionality depends only on 
the domain and the initial triangulation To- We define the restriction operator 
T^ih ■■ ^ V2h by 

{ihu,v)^ = {u,n2hv)2h > e V2h,v e Vh , 

where '■ V2h ^ V/i is the inclusion operator. It follows that 72.^^ is uniformly 
bounded with respect to /i G /, that is, there exists C independent oi h E I so that 

\\T^2hA^C\\u\\, VueVh. (4.14) 

We call a triangulation T locally symmetric if for every vertex P the union of 
triangles in T having P as a corner is invariant with respect to the reflection through 
P given by rp{x) = (2P — x). If T is locally symmetric and tpp is the nodal basis 
function at P, then tpp o rp = ipp. Furthermore, a simple calculation shows that for 
any linear map L{x) = aixi + 02X2 we have 

i Lpp{x)L{x~ P)dx^O . (4.15) 
Jo. 

Naturally, a uniform mesh is locally symmetric. 

The following grouped results are either simple consequences of Condition 12.11 or 
extracted from [T3j. 

Lemma 4.4. // £, Ch satisfy Condition \2.1\ there exist constants C{C) and 
C — C"(f2) independent of h such that the following hold: 

(a) H^,L^ - uniform stability of Ch: 

WChulH'^.^n) ^ C{C)\\u\\ ,yueVh,m = 0, 1; (4.16) 

(b) smoothing of negative-index norm: 

\\Cu\\ sC C{C) \\u\\ , Vu eVh,m^ 1,2 ; (4.17) 

(c) negative-index norm approximation of the identity by 7r2/i, 72.^^; 

!(/ - ^2/.)u||5_.(f,) s; c'h^ \\u\\ , yu e Vh; (4.18) 
\\{i-n^hMfi-.(n)^c'hP\\ui yueVh, (4.19) 
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where p — I on an unstructured grid, and p = 2 on a locally symmetric grid; 
(d) C diminishes high-frequencies: 

mi - n2h)u\\ ^ C{C)h^ \\u\\ , yueVh-, (4.20) 

\\c{i ~ n^Ml ^ cic)hp \\u\\ , yueVh , (4.21) 

where p — 1 on an unstructured grid, and p — 2 on a locally symmetric grid; 
(e) 

\{Cu,Cv) - {ChU,CHv)\ ^ C{C)h''\\u\\ ■ \\v\\ ,yu^Vh. (4.22) 

Proof. The stability conditions at (a) are direct consequences of (12.31) and (|2.4p . 
and (b) follows from (|2.3|) (see also 13 Corollary 6.2). The estimate ()4.18p is a 
straightforward consequence of the Bramble-Hilbert Lemma [9] , while (|4.19|) follows 
from Theorem 6.6 in [T^ (see Example 6.7 for the uniform mesh case). The inequalities 
at (d) follow from (b) and (c), and (e) follows from (|2.4p and the uniform boundedness 
of □ 

Proposition 4.5. If the operators C, Ch satisfy Condition 12.11 with the weights 
given by (j2.1l) . then Condition 14.11 holds with C{'H) ~ C{C) and p — 2 if the meshes 
are locally symmetric, or p — \ otherwise. 

Proof. To simplify notation we write tt = 7r2/i, TZ = 7?.^^. First we have 

{C2h^2h'^u,v)^ = {C2h^2hT^u,TZv)2h^ {^2h'^u,C2h'R-v)2h , and 

Therefore 

\{£*2h£2hTrU,v)i^- {ClChU,v)^\ ^ \{C2hTTU,C2h'R-v)2f^- {ChU, Chv) 
< \{C2hT^U,C2h'R'V)^^- {C2hT^U,C2h'R-v)\ + \{C2hTrU,C2h'R-v) - (ClTU, CTZv)\ 

^ y ^ 

Ai A2 

+ \{CT:u,Cnv) - {Cu,Cv)\ 

" v ' 

A3 

+ \{Cu,Cv) - {ChU,Chv) \ + \ {ChU,Chv) - {ChU,Chv)f^\ . 

^ ^ ^ V ^ 

A4 As 

For Ai we have 

l l4T3l igjej 

Ai ^ C{C)h^\\C2h7Tu\\Hi(n)-\\C2h'Jlv\\min) < CiC)h^\\Tru\\ ■ \\TZv\\ 
1I4T4I 

C{C)h^\\u\\ ■ \\v\\ , 
and a similar estimate holds for A5. Also (|4.22p implies that 

max(A2, A4) < C{C)h'^\\u\\ ■ \\v\\ . 

For A3 we have 

1 14:201 . don 

A3 ^\{C{Tr- I)u,Cnv)\ + \{Cu,{C{n- I)v))\ s; C{C)hP\\u\\-\\v\\ . 
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Since p ^ 2, ^3 is the weak link, and we have 

\{['H2h^-nh)u,v),} = \{{Cl,,C2h^-ClCh)u,v),} ^ C{C)h^u\\ ■ \\vl Vt; e V;, , 

and the conclusion follows from the equivalence of || • || with ||| • □ 
For p e W'^{Vl) we denote by 

\\p\w^(n)/R = max ||9"p||ioc(fi) , 

1^ |a| ^2 

which effectively is the norm on the quotient space W^(f2)/M. 

Lemma 4.6. If p & W^(f2) then there exists a constant C > independent of 
/9, h so that 

\\{Vp - V^p^)u\\ ^ Ch^\\p\\w-^i^)/Au\\H^{n) , e V,„ /i e / , (4.23) 
where ph = Ih{p)- 

Proof. Note first that V'^ = = Ih'Dp, since only the node values of p enter 
the definition of the operator V^. Given u GVh, for each triangle T S we have 



\pu\h^t) < 

because u is linear on T. Therefore 



W^{T)/R • |p||/i-i(T) ' 



(4.24) 



\{Vp-V'^)u\\ = \\pu-Ih{pu)\\^Ch^ I \pu\l2^T^ 



< Ch^ ( ^ IIpIIvk^(t)/bII 



^i(T) I ^Ch^\p\wu^)lv\Aii 



(n)/K||"l|ffi(f2)- 



Lemma 4.7. If p ^ W^{rt) then there exists a constant C > independent of 
p, h so that 

KPp -I?;,'Ju||5_,(^) < ChP\\p\\w^^n)/m\\u\\ , W&V,,, h&I , (4.25) 

where ph = 1-h{p), and p = 2 if the mesh is locally symmetric, otherwise p = I. 

Proof. We focus on the situation when the mesh is locally symmetric and leave 
the general case as an exercise. Let u € V/j,t; G IIq{^) be arbitrary, and denote by 
Si — supp((^^). The constant C is assumed to be independent of u,v,p,h. First 
remark that each triangle in Th lies in at most three of the sets Si and that 

diam(S'j) «C C/i, 1 i Nh (4.26) 
due to the quasi-uniformity for the meshes. Also note that 

p^'^-lHip^'^)=^'^ip-piP.h), l^^^N„. (4.27) 
Further we define Vi — ^^.^^^^-g^ v, for 1 ^ i ^ Nh, and Ui = u{Pl''). Then 
u = J2i^=i ^iVi and 



\{{V,~V'^Ju,v)\^ 



(pu -Ihifu)) V 

v^i{p~p{pmv-v.) 



4=1 



vi{p-p{ph)v 



S^ 



^^{p-p{Pl'))n. 
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For the first term in the sum above 



^ y'l\\-\\p-piPt)\\L^(S0-\\v-V,\\L2^S0 



< Ch'^hiW ■ |p|iyi(s.) ■ \v\h^s,) ■ 
For the second term in the sum we take advantage of the local grid symmetry 

h/ / T-,h\\ <4.15l 

'P^{p-p(Pn)v^ 



^ Ch'^bt\L^(s,)\W'l\\-\p\w^(s.) ■ 



Since ||tii||^2(5.-) ^ II''^IIl2(5.-) we now have 
(pu-Ihipu)) V 



< Ch'^WpWwim/mYl^Uil ■ y'lW ■ \\v\\hi(s,[ 



i=l 



(n)/R 



El 



/ill 2 



{S^ 



< Ch^\p\wi,{n)/Au\ ■ \v\H^{n) , 



where for the last inequality we used the quasi-uniformity of the mesh and the fact 
the each triagle is in at most three of the sets 5*1 . The conclusion follows after dividing 
by ||w||_ffi(n) and taking the sup over all v £ Hq{^). □ 

Proposition 4.8. If the operators K, Kh satisfy Condition 12.11 with v = 2 on the 
locally- symmetric meshes Th with the weights given by p. II) . then C, Ct also satisfy 
Condition 12.11 with 



C{C)^C{K.)\\\-^\\wu^) ■ 

If the meshes are not locally- symmetric, then the power of h in Condition 12.11 [b] for 
the operators C, Ch is 1 for both m = 0, 1. 

Proof Note that for p e L°°{n), \Dpu\ ^ ||p||L~(n) ' therefore 

||/C2?i/^u|U™(o) ^C(/C)|pi/^u|KC(/C)||A-^|U»(o)-||u|| , m = 0,l,2, 
which implies the smoothing condition (|2.3p . For u G 
||(/CI?,/^-/C,p;'/^)«||h.(o) 

< ||/C(2?i/vy - V\^^)u\\H^in) + ||(/C - lCn)V\^^^u\\H^n) 

EH T 1 1 

s$ C{]C)h'\\\ ^\\w2^(a)/AA\H^{a)+C{lC)h\\\-^\\L^(a)\\u\\ 

< C{K)h\\\-^\wu^)\\u\\ , 
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where in the last inequality we have used an inverse estimate. This proves (12.41) for 
/C,/C/i with m — I and concludes the proof for the non-locally-symmetric case. For 
m = 2 and locally-symmetric mesh 

^ Wm^i/vi - ^vta)^! + II(^ - 

BTft . drill 

sc: CilC)h^\\X-^\\wun)/M\\u\\ + CilC)h^\\X-^\\L^^n)\\u\\ 
^ C{IC)h'\\\-'^\\wun)\\u\\. U 

We conclude the analysis of the two-grid preconditioner with the computation 
of the spectral distance (SD) between {Gh)~^ and the two-grid preconditioner Sh 
defined in (|4.8p . This step facilitates a smooth transition from the two-grid analysis to 
the multigrid analysis of the next section. Given a Hilbert space {X, (•,•)) we denote 
by the set of operators with positive definite symmetric part: 

£+ (X) = {T e £{X) : {Tu, u) > 0, Vu G A" \ {0}} . 

First we define the joined numerical range oi S,T £ £^{X) by 



[ {Tew, w) J 



W{S. 

where Tc{u + iv) — T{u) + iT(w) is the complexification of T. Note that if T is 
symmetric positive definite, then W{S, T) is simply the numerical range of T-^ST-*. 
The spectral distance between S,T £ introduced in ^13j as a measure of 

spectral equivalence between S and T, is defined by 

d^(S',T) = sup{ I In z I : zeW{S,T)}, 

where In is the branch of the logarithm corresponding to C \ (— oo,0]. Following 
Lemma 3.2 in [ig, if W(5,r) C 6„(1) = {z e C : |z - 1| < a} with a e (0, 1), then 

d.{S,T) sup{|z-l| : zeW{S,T)}, (4.28) 

a 

which offers a practical way to estimate the spectral distance when it is small. The 
spectral distance serves both as a means to quantify the quality of a preconditioner 
and also as a convenient analysis tool for multigrid algorithms. Essentially, if two 
operators S, T satisfy 



1 - 5 



{Sew, 



{Tew, w) 



s$ l + (5 , Vwe A'^\{0}, 



with 6-^1, then dcr{S,T) Ri 6. li N Ri G ^ is a preconditioner for G, then both 
da{N,G~^) and daiN~'^,G) (quantities which are are equal if G, N are symmetric) 
are shown to control the spectral radius p{I — NG) (see Lemma fA. 2 1 in Appendix lAl for 
a precise formulation), which is an accepted quality-measure for a preconditioner. The 
advantage of using d„ over p{I — NG) is that the former is a true distance function. 
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Theorem 4.9. // the operators K,, tCh satisfy Condition 12.11 on the locally- 
symmetric meshes Th with the weights given by (12. ip . and A G W^{fl) satisfies 
Z/iA = Xh, there exist C,S > independent of h, A so that for h'^\\X~'^\\w^{n) ^ ^ 

{GJ;\Sh) < Ch^\\X-^^\\wl(n) , (4.29) 

where Qh and Sh are defined as in (|4.4|) and ()4.8|) . // the meshes are not locally- 
symmetric then the power of h in (|4.29p is 1 . 

Proof. Again, we restrict our attention to the locally-symmetric case. The oper- 
ator Qh is symmetric positive definite with respect (•, •)^ and satisfies 

{Qhu^u)^^ Ml + \\Chuil^ Ml . 

Therefore <j{Qh) C [l,oo), and 1 for all < 0. Due to the norm equivalence 

III -I ft I'll there exists Ci > so that ||CJ^ ^ || ^ Ci. By LemmalT^and PropositionsH31 
and 14.81 we have 

\\gh-Mh\\^C2h^\\\-Hwun) 
for some constant C2 > 0. Hence for C3 = C1C2 

\\i - g-JUhghh\ ^ WQ^V ■ WQh-m ^ C3/i2||A-5||^^(f,) . 

Since for any operator T we have ||/ - T^^l (1 - a)-^\\I - T|| if ||/ - T|| ^ a, 

1 1 Sh=J^^^ 4 _i _i 4(7, ^ 1 
WI-QIShQIW ^ ^\I^Qn''Mhgn~'H^h^\\\-n\wu^) , (4.30) 

provided ft,^||A^^ || ^2 (j^) ^ (5 = By further restricting 5, we can assume the right- 
hand side of (1430)) to be 1/2, which implies that WiSh^G^^) C Bi{l). By (1428)) 
we obtain 

d,{Su,g^^) ^ ^C3/i2||A-i||H^^(n) . □ 

The highlight of the last result is the presence of 0{h^) (or 0{h) for general 
quasi-uniform meshes) in the right-hand side of ()4.29p . which is the optimal order 
of approximation in h. We should stress that for classical multigrid methods for 
differential equations one has 0(1) as the right-hand side estimate, which is sufficient 
for mesh-independence. In this case, if the theoretically introduced smooth function 
A could be the same for all meshes, the number of iS/j-preconditioned iterations is 
expected to decrease with h I 0. In reality, the discrete function A^ is tied to the 
Lagrange multipliers vi,V2, which in turn are related (actually expected to converge 
to as ^,h l 0) the Lagrange multipliers A, A of the continuous problem. Since in 
general the latter are only in L^, the factor ||A~ 2 [^a^^f^^ is expected to be unbounded 
as ^ 4, 0. Therefore the preconditioning qualities of Sh are expected to increase 
with h I 0, but decrease with | 0. Thus for large-scale, high-resolution problems, 
where h <^ 1, the presented method is expected to perform very well, especially in 
connection with the multigrid method discussed in the next section. However, for 
fixed ft., as /i J, in the IPM formulation and the approximate solution approaches u, 
if the inequality constraints are active then the quality of the proposed preconditioner 
normally degrades. The advantages or disadvantages of this method will ultimately 
be discussed based on numerical experiments in Section [Sj 
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5. The multigrid preconditioner. While the two-grid preconditioner Sh may 
be efhcient in terms of number of iterations, it is expensive to apply. In this section 
we develop a multigrid preconditioner Sj^^ that also satisfies the optimal order es- 
timate (|4.29p but has a lesser cost. Since the process of passing from a two-grid to 
a multigrid preconditioner of comparable quality has been analyzed in [13], we give 
here only a brief description. In this section we assume a finite number of grids 



— z 



and the goal is to ultimately construct an efficient multigrid preconditioner for the 
operator on the finest grid. 

Consider the operator : £,{Vhi_i) — > ■C(Vh.) by 

Cf. (|4.8p we have Sh^ = Xi_i{Gf^l^_^). If we define Sj^, h € /max, recursively by 

then has a V-cycle structure. However, it is shown in jl3] that is suboptimal, 
in that it satisfies (|4.29l) with h'^ replaced by /ig. Thus the quality of does not 
improve with ft, | 0, as desired, it is simply mesh-independent (it only depends on ho). 
To achieve the desired result we define the operator 9Ti : £(Vh.) £(V/iJ 

The latter is related to Newton's method for the operator-equation — Qh^ = 0; 
namely, if Xq is a good guess at the solution, i.e., approximates well 5^.^, then the 
first Newton iterate starting at Xq is Xi = ^i{Xo) (see also Remark 3.11 in [IS])- We 
define the multigrid preconditioner using the following algorithm: 

Algorithm 1: Operator-form definition of S^.^ 



2. 5™" := Qj^ 7. coarsest level 



1. if i = 

3. else if i < imax 

4. '^i{3\_^{Sl\^^J) 7. intermediate level 

5. else 

6. at 1 (5™^^ ) 7. finest level 

7. end if 

8. end if 

The key factor in Algorithm 1 is the application of 'Tli at Step 4, and here is why: 
while Q^^ is well approximated by 3-_;^(5^"^J provided that Ghi_^ ~ ^h^-i (recall 
that Q~^KiS]n = 3(S/7i^i))' apphcation of Vli brings 3(5^^^) even closer to 
This step is critical if we want dcr{G^^ ,Sj^^) — 0{dcr{G^^ ,Sh))- Also, there are two 
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main reasons for splitting the cases of intermediate vs. finest resolution, as opposed 
to just replacing with J) in (|5.ip . First we would like to 

have = Sh for h = /ii,-„ax if only two grids are used. Second, the application of 
Uli includes a multiplication by Ghi since for the intended large-scale applications the 
finest-level mat-vec is expected to be very costly, we prefer that no such mat-vecs are 
computed inside the preconditioner. 

Theorem 5.1. In the hypotheses of Theorem 14.91 and with 5™^ defined as in 
Algorithm 1 there exist C,S > independent of h and \ so that for h'^\\X^^\\]Y2 (qj ^ S 

da {G^\S["') < Ch^\\X-i\\wiin) , ioTh = /i^ax . (5.2) 



The proof of Theorem 15.11 follows closely that of Theorem 5.4 in [13] and, in the 
interest of brevity, we do not give further details. SufBce it to say that the use of the 
spectral distance is instrumental, and that an essential ingredient is the symmetry 
(with respect to (•, •)^) oi Gh- 

In practice, for large-scale problems, neither Gh nor 5™^ are ever formed, so 
both are to be applied matrix-free. A simple verification shows that, given some 
M e £(Vft.), the vector ti = {^i{M))r can be computed by setting u U2 where 
Uk+i '■— Uk + Mir — GhiUk) with uq = 0. Thus the matrix-free application of is 
computed by the following function; 

Algorithm 2: Matrix-free implementation of the action u — S'^^r. 

1. f unction M = AfG'(r, z) 

2. if 1 = °L coarsest level 

3. u :— Gh^r °L direct or unpreconditioned CG solve 

4. else 

5. u:={I- nh^_,)r + MG{nh^_,r, i - 1) 

6. if i < iaia.K 7o intermediate level 

7. ri := r - GhiU 

8. Ml ;= (/ - TThi^Jri + MG{Trh,_^ri,i - 1) 

9. u := u + Ul 

10. end if 

11. end if 



As can be readily seen. Algorithm 2 has a W-cycle structure. To estimate the cost 
of MG(-, imax) we denote by T{i) the cost of applying MG{-,i) for < z < imax- If 
we assume that one residual computation at level i has complexity 0{Nf^\ and that 
the cost of computing an L^-projection is negligible compared to that of a residual 
computation (this is reasonable for most applications since mass matrices are normally 
easy to invert) then the resulting recursion for the function T reads: 



m = 0{Nl) + 2T[i-l) . 

For i = i„iax the term 0{N'^,) is replaced by a potentially smaller cost of just comput- 
ing an L'^-projection. Given that Nhi_i ~ 2~'^Nhi (here d = 2), a standard argument 
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shows that 



Tit^,^)=0{NlJ , 



that is, a cost that is proportional to that of a residual computation. 

Another comment refers to a detail that is not very transparent in Algorithm 2, 
namely that coarse versions of Xh^ are necessary for each level, because 



Since the original problem is solved starting at the finest level, where Xh^^^^ is given 
by the optimization algorithm, the functions Xhi are obtained by simply discarding 
the values at finer nodes of A/j^^^^. The parameter /? is hidden in A and affects the 
process indirectly. 

6. Applications and numerical examples. In this section we discuss two 
applications. The first is related to the inverse contamination problem studied in [2j[T], 
where /C is a time-T solution operator of a parabolic equation. The second is a 
standard elliptic-constrained optimal control problem with additional box-constraints 
on the control. 

6.1. Solution strategies and metrics for success. For both applications we 
apply Mehrotra's algorithm and we solve the inner linear systems using the multigrid 
preconditioner previously defined. In the absence of multiple grids, the linear sys- 
tems (|3.12p are solved using conjugate gradient (CG), while for more than one level 
we used MG-preconditioned conjugate gradient squared (GGS), because of the slight 
non-symmetry of the MG preconditioner. As a first metric we record the number of 
inner linear iterations needed at each outer iteration; secondly, we record the total 
number of finest-grid mat-vecs for the entire solution process. Recall that each outer 
iteration requires the solution of two linear systems with identical matrices, namely 
one for the predictor step and one for the corrector step; in the interest of the pre- 
sentation we record only the linear iterations for the predictor step. Also, a small 
number of mat-vecs are required in the process in addition to those needed for the 
predictor-corrector solves, and are refiected in the count. With respect to the second 
metric we remark that the proposed algorithm is intended for large-scale problems, 
with the most expensive computation being the finest-scale residual computation. 
The ultimate goal is to significantly reduce the total number of finest-scale mat-vecs, 
because this is expected to be directly linked with execution time in a truly large- 
scale computation. With regard to the second metric we would like the total number 
of finest- level mat-vecs to decrease with h ], 0. As for the first metric (number of 
iterations) we would like to witness the following: 

[a.] The number of MG-preconditioned CGS-iterations should be less than half of the 
unpreconditioned CG-iterations (each CGS-iteration involves two mat-vecs, while a 
CG iteration requires only one). 

[b.] For a given resolution h, the number of MG-preconditioned GGS iterations should 
be relatively bounded with respect to the number of levels used, provided the coarsest 
level is sufficiently fine, as stated in Theorem 15. II 

[c] Mostly important, the number of MG-preconditioned CGS iterations 
should decrease with /i ^ 0; in other words, the MG-preconditioned GGS becomes 
increasingly advantageous compared to GG as the problem-size increases. One word 
of caution though: linear systems of different resolutions are not necessarily related in 
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a direct fashion, since their "A" is dictated by the evolution of Mehrotra's algorithm, 
which is slightly different for each resolution. For example, the tenth linear system to 
be solved in the IPM process at a resolution 2h is not necessarily some coarse version 
of the tenth system to be solved at resolution h. 

Also, we should point out that in all our tests we use a cold start, that is, we 
do not take advantage of results from coarser levels except for in the MG-solve of 
the inner linear systems. While for realistic applications "warm-start" strategies are 
essential, we restrict our attention to the way our multigrid technology plays a role 
in solving the inner linear systems. 

6.2. Time-reversal for a parabolic equation. We consider the problem of 
finding the initial state (the control) for a system governed by a parabolic equation 
given the state at a later time T under additional box-constraints on the control. 
Multigrid-preconditioning for the unconstrained version of this problem was studied 
in detail in [TH [13] and for the space-time measurements in [5j . 

Formally, we consider the following parabolic initial value problem with periodic 
boundary conditions 

r dty - d,{ad,y + by) + cy = , on [0, 1] x (0, T] , 

i y{0,t) = yil,t), dMO,t) = d,yil,t) , foriG(0,r], (6.1) 

[ y{x, 0) — u{x) , for a; G [0, 1] , 

where a > 0, 6, c are constants, and T > is the end-time. For i > we denote by 
S{t) E £(L^([0, 1])) the time-i solution operator mapping the initial value onto y{-,t) 

u I — > y(-,t} , 

and let /C = S(T). The discrete ICh is obtained by using a Galerkin formulation with 
continuous piecewise elements on a uniform grid for the spatial discretization and 
Crank-Nicolson in time. It is shown in [30] (see also [36]) that for sufficiently small h 
the following estimate holds: 

\\ICu - IChTThulH'^ < Ch^-"'\\u\\ , Vu e L^{[Q, 1]) , m = 0, 1, (6.2) 

where C — C'(T), provided the time step k is proportional to the spatial resolution 
k — Cih, with C'l chosen to ensure stability. Consequently, space and time resolutions 
are refined at the same rate, and Condition 12. II is verified, so our theory applies. The 
specific details (boundary conditions, constant advection etc.) in this example were 
chosen for convenience, however, two and three spatial dimensions, other types of 
boundary conditions, as well as smoothly varying functions in place of the constants 
a, 6, c are supported. 

A direct verification of the convergence order. As mentioned earlier, when running 
Mehrotra's algorithm with different resolutions, the added diagonal terms A may not 
be in direct relationship with each other. Hence, in order to practically verify the 
presence of in the estimate (|4.29|) we resort to an artificial context: we construct 
Gh based on a fixed function A/i ~ I/i(sin) -I- f3 for hj = 80 • 2^,j = 0, 1,2,3, and 
we define the corresponding two- grid preconditioners Mh- Then we compute the 
"distances" dh = max{|lnQ;| : a G iy{Ghi-^h)}- Since dh approximates the spectral 
distance of interest (because A/^ is close to being symmetric, actually we have dh ^ do-) 
we expect to see that dh — 0{h?). We repeat the experiment for /3 = 1, 0.1, 0.01. The 
results presented in Table [131 while not yet converged, give a strong indication of an 
asymptotic rate lim/i_j.o d2h/dh = 4. 



24 



A. DRAGANESCU AND C. PETRA 



Table 6.1 

dfe = max{|lnQ!| : a £ u{Qfi,Afii)} for \{x) (sin(a;) 



h\l3 


1 


0.1 


0.01 




dh 


rate 


dh 


rate 


dh 


rate 


1/80 


0.0206 




0.1127 




0.2812 




1/160 


0.0066 


3.1342 


0.0363 


3.1078 


0.1270 


2.2140 


1/320 


0.0020 


3.3140 


0.0102 


3.5488 


0.0445 


2.8535 


1/640 


0.0006 


3.5199 


0.0027 


3.7365 


0.0123 


3.6284 



Numerical study. We consider the "true" initial value uq supported on two in- 
tervals with Mo reaching the value 1 on one of the intervals and 1/2 on the other 
interval, then let / = ICuq. Specific values are a = 4 • 10~^,& = 0.4, c = 0, and 
T = 0.8. In Figure [Ql we show uq,/ as well as the converged solution Mmin of the 
box-constrained optimization problem with /3 = 10"'^, a value chosen because of the 
relatively good (visual) agreement of uq with Umin- We run Mehrotra's algorithm for 
h = 2^^", 2^^^, 2^^^, 2^^^, and for each h we test the solvers with 1, 2, and 3 levels. 
The number of linear iterations required by each of the linear solves in the predictor 
step are shown in Figure [6?2| while the corresponding values for ||A~^/^||vi/2 and /i are 
shown in the top two pictures of Figure 16.31 




Fig. 6.1. Solution with f = Kuq, /3 = 10~^, s£ u 1. 

First we remark that the number of unpreconditioned CG iterations appears to 
be mesh- independent (top chart in Figure 16. 2p : essentially the curves representing 
the number of iterations for each of the resolutions more-or-less overlap. We notice 
only a slight increase in number of iterations for higher resolutions. Second, from the 
middle chart in Figure 16.21 representing the number of two-grid CGS iterations, we 
infer that the number of two-grid preconditioned iterations consistently decreases 
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Fig. 6.2. Number of iterations for each of the predictor-step linear systems solved (/3 = 10 
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Fig. 6.3. Top: discrete ||A^ ^ based on discrete Laplacian should give an idea of the size 

of ||A~2 jjj^2 ; middle: log-plot of fi as a function of outer iteration; bottom: X/^ for h = 1/4096 at 
the last outer iterate (/3 = 10~^). 



with h 0, as desired. For example, at the twenty- fourth iteration these numbers 
are 8,6, and 4, while at the twentieth they are 7,4,4. This phenomenon is repeated 
for the three-grid preconditioner as can be seen from the bottom chart in Figure 16.21 
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Moreover, for the most part, the number of MG-CGS preconditioned iterations is 
significantly smaller than half the number of unpreconditioned CG iterations: e.g., for 
h = 1/4096, up to the thirteenth iteration (where fi is already down to approx. 10~^, 
see Figure |673|) only one two-grid (or three-grid) MG-preconditioned CGS iteration is 
necessary, while up to 15 unpreconditioned CG iterations are needed. However, after 
the thirteenth iteration ||A~^/^||w^ shows a significant increase, as seen on the top 
chart of Figure 16.31 and the MG-preconditioned CGS is less efficient: for h = 1/4096, 
at the twenty-fourth iteration 4 two-grid iterations are needed compared to 13 CG 
iterations, a lesser advantage compared to the earlier outer iterations. The bottom 
chart in Figure 16.31 shows the last computed A^, at /i = 1/ 4096 to give an idea of why 
the quantity ||A~^/^||vv'2 is so large. A comparison between the bottom and middle 
charts in Figure 16.21 shows that the number of MG-preconditioned iterations is not 
very sensitive to the number of levels, provided the coarsest mesh is sufficiently fine. 
In this example four levels would force a much too coarse base mesh, and produce 
unsatisfactory preconditioners. The last piece of evidence is the total count of finest- 
level mat-vecs, shown in Table 16.21 In this example, a mat-vec involves solving the 
advection-reaction-diff'usion equation on [0, T]. The data clearly shows that, as h 10, 
the two-level solvers is getting increasingly efficient in this metric compared to CG: 
the ratio goes from 581/728 for h = 1/1024 to 377/768 for h = 1/8192 = 2-^^. We 
should remark also that the essential impediment to a more significant improvement 
over CG lies in the increase in the ||A~^/^||n/2 as h ^0. As shown in Theorem 14. 9i 
non-smoothness of A~^/^ decreases the preconditioner's eflaciency. 



Table 6.2 Table 6.3 

Total number of fine- grid mat-vecs Total number of fine- grid mat-vecs for 

for the ID reversed parabolic equation the 2D elliptic- constrained opt. ctrl. problem 



h \ levels 


1 


2 


3 




h \ levels 


1 


2 


3 


4 


1/1024 


728 


581 


661 




1/256 


354 


282 


572 




1/2048 


740 


463 


489 




1/512 


355 


220 


250 


452 


1/4096 


764 


403 


425 




1/1024 


355 


198 


210 


224 


1/8192 


768 


377 


403 




1/2048 


363 


172 


174 


174 



6.3. An elliptic-constrained control problem. In this example we discuss 
the elliptic-constrained optimal control problem (jl.4|) from Example B, which is a stan- 
dard test problem in PDE-constrained optimization ^ [26] , and corresponds to (|l.ip 
with K- — . We consider a square domain with a continuous piecewise linear finite 
element discretization based on the standard three-lines triangular mes Standard 
estimates for finite element solutions of elliptic problems show that Condition 12.11 is 
verified |9]. 

Numerical study. Let 51 = [0,1] x [0,1], /? — 10^^, and / be the function that 
satisfies A/ = uq, /|ao = 0, where uo{x, y) = | sin(27ra;) sin(27ry). With this selection 
of /, the choice u = uq would be a solution of (|1.4I) if /3 = and no box constraints were 
present (or if [— 1,|] C [m,u]). Here the bounds [u,u] = [—1,1] are active: without 
them, given that /3 <C 1, the solution would be close to uq, that is, would have a 
maximum (resp. minimum) close to 3/2 (resp. —3/2). The solution with h = 1/128 
is depicted in Figure We have solved the problem with h = 2^^, 2^^, 2^^°, 2^^^ 



^The three-line mesh is obtained by dividing the square into equally sized squares with sides 
parallel to the coordinate axes, and by further cutting each little square along its slope-one diagonal. 
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using one, two, three, and four levels (where appropriate) using the strategy described 
in Section 




Fig. 6.4. Solution with f satisfying Af = | sm(27rx) sin{2ny), jS = 10^^, [u,u] = [-1,1]. 
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Fig. 6.5. Number of iterations for each of the predictor- step linear systems solved {p = 10 ^). 



As with the previous example, we show in Figure 16.51 the number of iterations 
required by each of the linear systems at the predictor step. The top plot shows 
the number of unpreconditioned iterations to level off at 21. The middle plot again 
shows two facts: the number of MG-preconditioned CGS iterations decreases with 
h \. Q. In addition, for the finest grid, the number of two-grid preconditioned CGS 
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iterations is less than 1/4 that of the number of iterations for the unpreconditioned 
case even when looking beyond the eighth outer iteration, where as before, roughness 
of the A-function lowers the quality of the MG-algorithm. For example, for the tenth 
outer iteration with h = 1/1024, CGS required 3 iterations, while 21 iterations were 
needed in the unpreconditioned case; for the eleventh outer iteration the numbers are 
5 vs. 20. Of course, each two-grid preconditioned CGS iteration is significantly more 
expensive than an unpreconditioned CG iteration, which is why the number of levels 
should be maximized. The bottom plot in Figure 16.51 shows that the two-level be- 
havior is replicated using three-level preconditioners. Moreover, at fine resolutions 
{h ^ 2~^°), the numbers of required three-level preconditioned CGS iterations are 
not significantly higher than those of two-level CGS iterations. However, this is not 
the case for low resolution h = 2~®, where the number of three- level CGS iterations 
required for the last three systems (not shown on the plot) is quite large: 45, 17, 58. 
This is why we insist that the MG-preconditioner is efficient only when the coarsest 
resolution used is sufficiently fine and the finest resolution h is small. With respect 
to the second metric, we show in Table [6731 the total number of fine mat-vecs for each 
of the runs, and the results confirm that the MG-preconditioner becomes increasingly 
efficient with h I 0. Recall that for this application a mat-vec requires solving the 
Poisson equation. 

Appendix A. Some facts about the spectral distance. Throughout this 
section (A", (•,•)) is a real, finite dimensional Hilbert space with norm ||- 1|. All operators 
in this section are assumed to be in C+ {X) (see Section 14.21 for definition) unless 
otherwise specified. The following inequalities were proved in [13) (Lemma 3.2): 

Lemma A.l. If a e (0, 1) and z e Bail), then 

ln(l -|- a) , , ,, , |ln(l — q:)|, , ^ 

^ Inz^^ ^1-^ • (A.l) 

a a 

For |lnz| ^ S we have 

i^^llnzK |l-zK ^^|lnz|. (A.2) 





Lemma A. 2. Let L,G e £+(<%") such that 

min(d,(L-\G),d,(L,G-i)) ^ (5 . 

Then 

p{I - LG) ^ ^—mm{d,{L-\G),d,iL,G-^)) . (A.3) 



Proof. If A G o'(/ — LG) then there exists a unit vector u £ X'^ such that 
(/ — LG)u = Xu, therefore 

(1 - \)u = LGu . (A.4) 
After left- multiplying with and taking the inner product with u we obtain 

(1 - A) (L-^u, u) = (Gu, u) , therefore A = 1 - /'^"''"^ . 
^ ' (L^'-u,u) 
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If we substitute v — G in (|A.4[) and take the inner product with v we have 

(1 - X)G-^v = Lv . therefore A = 1 - . 

(G ^v,v) 

Hence, if d„{L^^,G) ^ 6, then 

p{I-LG) sC sup{|l-z| : z = (Gw,u) /<L"^u,m) for some u e A-^XIO}} 

EI} - I 

< '-—d„{L-\G). 



Instead, if dcr{L,G^^) ^ 6, then 

p{I~LG) sC sup{|l-z| : z= (Lit, m)/(G-^u,w) for some ueA'^\{0}} 

< —^da{L,G~') . 
which proves (|A.3p . □ 
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